Asymmetric velocity correlations in shearing media 
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A model of soft frictionless disks in two dimensions at zero temperature is simulated with a 
shearing dynamics to study various kinds of asymmetries in sheared systems. We examine both single 
particle properties, the spatial velocity correlation function, and a correlation function designed 
to separate clockwise and counter-clockwise rotational fields from one another. Among the rich 
and interesting behaviors we find that the velocity correlation along the two different diagonals 
corresponding to compression and dilation, respectively, are almost identical and, furthermore, that 
a feature in one of the correlation functions is directly related to irreversible plastic events. 



I. INTRODUCTION 

In collections of particles with repulsive contact inter- 
action there is a transition from a liquid to an amorphous 
solid state as the volume fraction increases — the jam- 
ming transition. It has been suggested that this transi- 
tion is a critical phenomenon with universal critical expo- 
nents [1] and the successful scaling of rheology data from 
simulations is strong evidence that that actually is the 
case [2-4]. The precise values of the critical exponents, 
however, continue to be a matter of discussion [5]. 

At the very heart of critical phenomena is the notion 
of a correlation length that diverges as the critical point 
is approached. It is therefore important to identify the 
proper correlation length. Several works have tried to 
look for a growing order in the static quantities, but 
without much success. Another possibility is to look for 
a growing length in the dynamics. Velocity correlations 
in sheared systems were studied in Ref. [6] though not 
revealing any growing length. A large correlation length 
was however found in Ref. [7] and they also argued for a 
pronounced angular dependence of the velocity correla- 
tions [8]. 

In a previous work we reported the finding of a grow- 
ing characteristic length from the transverse component 
of the velocity correlation function [2]. The extraction 
of the correlation length exponent, however, seems to be 
more complicated than presented there and we therefore 
set out to do a more thorough analysis of the velocity 
correlations. As an important step in that direction we 
here consider some symmetry properties of velocity cor- 
relations in a sheared system and find a surprisingly rich 
and interesting behavior 

When shearing simulations are done slowly enough it 
becomes possible to separate the time evolution into elas- 
tic parts where the energy slowly increases and plastic 
contributions which are irreversible processes where the 
system rapidly evolves and dissipates energy [9] . We will 
argue below that the contribution from the plastic pro- 
cesses also may be seen in the velocity correlation func- 
tion. 

The content of the present paper is the following: In 
Sec. II we briefly describe the model and the simulations. 
Sec. Ill describes some rather direct measures of veloc- 



ity correlations and how they depend on the direction 
of the separation between the particles whereas Sec. IV 
deals with a more involved correlation function designed 
to capture the difference between clockwise and counter- 
clockwise rotations of the velocity field. A summary and 
some concluding remarks are given in Sec. V. 



II. MODEL, SIMULATIONS, AND MEASURED 
QUANTITIES 

A. Shearing dynamics 

Following O'Hern et al. [10] we simulate frictionless 
soft disks in two dimensions using a bi-dispersive mix- 
ture with equal numbers of disks with two different radii 
of ratio 1.4. Length is measured in units of the small 
particles (dg = !)• With rij for the distance between the 
centers of two particles and dij the sum of their radii, the 
interaction between overlapping particles is 
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We use Lees-Edwards boundary conditions [11] to intro- 
duce a time-dependent shear strain ^ — t'j. With peri- 
odic boundary conditions on the coordinates Xi and yi 
in an L X L system, the position of particle i in a box 
with strain 7 is defined as r^ = {xi + jyi, yi) which thus 
gives a shear flow in the x direction. We simulate over- 
damped dynamics at zero temperature with the equation 
of motion [12], 
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which is integrated with the second order Heuns' method. 

This above expression is for the total velocity, includ- 
ing the shearing part. In the analyzes of the velocity 
correlations below we will use the non-affine part of the 
velocity excluding the trivial shearing part yijx. This 
non-affine part of the velocity will be denoted by v. 

Our simulations are performed with N = 65536 with 
shear rates down to 7 = 10^*. The averages are typically 
from simulations during 1-3 days and nights with 128 
cores on a massively parallel computer. 
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FIG. 1. Configuration witli 4096 particles color coded accord- 
ing to the elastic energy of each particle. The dark particles 
(high elastic energy) make up force chains that preferably ex- 
tend along the x — y direction. 



III. VELOCITY CORRELATIONS 
A. Symmetry of a shearing system 

Figure 1 shows the presence of force chains in our sys- 
tem. The figure is a configuration with 4096 particles, 
color coded according to the elastic energy of each par- 
ticle. Note that the force chains [8] tend to be along the 
X — y direction, which is the direction of compression. 
This means that the force chains break the reflection 
symmetry along x and that the system is only symmetric 
under the combined transformation x — >■ —x and y — >■ —y. 
The same conclusion is readily drawn from the expression 
for the shear stress. 
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FIG. 2. Measure of the anisotropy in both particle velocities 
and contact forces from simulations with 7 — 10~^. Panel 
(a) shows the normalized (vxVy) for individual particles ver- 
sus density. The correlation changes from negative at low 
densities to positive around <j) ^ <j)j (dashed vertical line), 
signifying a change from a predominance of motion along the 
force chains to a slight overweight for motion perpendicular 
to the force chains. Panel (b) shows that the corresponding 
contact force correlation is always negative, which is consis- 
tence with the existence of force chains at all densities. This 
contact force correlation also has a marked feature aX 6 ^ Sj. 



where L is the linear system size, the sum is over pairs 
of particles, ff: is the x-component of the force between 
particles i and j, and yij is the y-component of their sep- 
aration. Note that the shear stress changes sign under 
the transformation x ^ —x but remains unchanged un- 
der the combined transformation x — > —x and y — > —y. 



B. Single particle properties 

As we will see below several symmetries that hold in 
systems at equilibrium are broken in a shearing system. 
The simplest symmetry is however respected; there is 
no net velocity in the system, (v) = 0. As remarked 
above, v is the non-affine part of the velocity. Here 
and in the following (...) represents the average over 
all particles and a large number of configurations gen- 
erated with our shearing dynamics. For the average ve- 
locity things are unusually simple since the same result 
holds for each individual configuration as a consequence 



of the overdamped dynamics and total force balance, 
V = E. Vi = CY^ij Uj = as fy = -fji. 

The conclusion of a vanishing average velocity may also 
be reached from the symmetry considerations. Since the 
combined transformation also implies the change of sign 
of both velocity components, Vx — > —Vx and Vy — )■ — Wy, 
it follows that (u^) — (— w^), for fi = x,y, which gives 
M = 0. 

In contrast to equilibrium results from symmetry that 
(vxVy) — 0, one finds that this quantity does not vanish in 
the sheared system. Fig. 2 shows (vxVy) / (v^) against (j). 
At low densities the correlation is negative which means 
that the particles tend to move more along than perpen- 
dicular to the force chains. The correlation changes sign 
at « 0.81, reaches a peak at (/> « 0.84 « (f)j and then 
decrease towards zero. This means that there is a region 
around (jjj where the particles move slightly more in the 
direction perpendicular to the force chains. 

A related quantity is {ffjlL) which is the correlation 
between the different components of the contact force iij . 
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FIG. 3. The anisotropy measured through the (vy) relative 

to (v^) — (vV) + {%)■ This fraction is close to 50% in the 
regions around jamming but decreases slightly with increasing 
density. 



FIG. 4. Total velocity correlation along the four different 
main directions. The decay of the correlations is monotonic 
along the main directions, x and y and non-monotonic along 
the diagonals. 



Note that there is a direct relation between the velocities 
and the contact forces: v^ = C^ ■ %, where the sum ex- 
tends over all particles j in contact with i. Nevertheless, 
{f'iifii) behaves rather differently from the velocity cor- 
relations. Whereas (v^Vy) / (v^) has a peak at (/) ~ (j)j, 
Fig. 2(b) shows that the normalized {ffjfij) changes at 
the same density from a rapid increase (which is a de- 
crease in magnitude) to being almost constant. 

Another measure of the asymmetry is the relative mag- 
nitude of the two velocity components. From energy 
balance — that the dissipated power has to be equal to 
the supplied power — follows the result for the velocity 
squared: C(v^) = {L^/N)aj [6]. This dissipated power 
needs however not be equal in the x and y directions. 
Fig. 3 shows that the fraction of the power dissipated by 
velocities along the y direction is close to 50%, but also 
that there is a clear dependence on density: (w^) / (v^) 
decreases from 0.506 to 0.495 when the density increases 
from 6 = 0.80 to 0.94. 



C. Spatial dependence 

We now turn to the spatial velocity correlations, i.e. the 
correlations between pairs of particles with separation r. 
We first focus on the total correlation function. 



5(r)^(v(r')-v(r' + r)) 



(1) 



and examine how it depends on both magnitude r = 
|r| and direction of r. We will study this correlation 
function in four different directions: along the two main 
directions, x, y, and along the diagonals 

s=ix + y)/V2, 
i={y- x)/V2. 

In our shearing geometry t is the direction of uniax- 
ial compression and s the direction of uniaxial dilation. 



Figure 4 shows the correlation functions along these four 
different directions. The curves are pairwise equal with, 
on the one hand, the correlations along the main direc- 
tions, x and y, and on the other hand the correlations 
along the diagonals. These results should be essentially 
without finite size effects since the system size is L w 300 
whereas r only extends up to 50. In view of the density 
dependence in Figs. 2 and 3 we have confirmed that the 
general behavior remains the same for a wide range of 
densities around (j)j. 

The result that g{r) behaves the same along both di- 
agonals is different from the earlier finding of an angular 
dependence of the correlation length in shearing systems 
[7] with minimum and maximum along the different diag- 
onals. The reason for this difference is not clear, but we 
speculate that it is related to the very different dynamics 
in their system, which is also reflected in the oscillatory 
behavior of their velocity correlation function. 

To investigate the reason for the dependence of g(r) on 
the direction of r, we separate the correlations into longi- 
tudinal and transverse components, parallel and perpen- 
dicular to the separation, respectively. In the x direction 
these components are 

9\\{rx) = {vx{0)vx{rx)) , 
9±irx) = {vy{0)vy{rx)) , 

with obvious generalizations to the other directions. Af- 
ter this separation we find that the transverse compo- 
nent along these four different directions behave about 
the same, see Fig. 5(a). (The rather small differences 
in the transverse component for the two diagonal direc- 
tions will be discussed further below.) The difference is 
largely due to the decay of the longitudinal component, 
see Fig. 5(b) which is monotonic along the main (x and 
y) directions but non-monotonic along the diagonals. 

The result above is a different behavior along the di- 
agonals compared to the main directions. We now in- 
stead focus on the difference between the two diagonal 
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FIG. 5. Transverse and longitudinal velocity correlations, re- 
spectively, along the four different directions. The transverse 
correlation in panel (a) is very similar for all four directions 
whereas the longitudinal correlation function shown in panel 
(b) is different for the diagonal directions, compared to the 
two main directions. The differences in the total velocity cor- 
relations in Fig. 4 are thus related to differences in the longi- 
tudinal correlations. 



directions, s and t. From Fig. 4 we found that the to- 
tal velocity correlations along these two directions are 
very similar. Nevertheless, as shown in Fig. 5 both 
g\\ and 5^ are clearly different. This suggests that 
this difference originates from the mixed correlations 
gxy{i^) — {vx{0)vy{'r)) ■ To See this we use the defini- 
tions g\\{rs) = (ws(O)us(rs)), and g±{rs) = {vt{0)vt{rs)) 
together with Vg — {vx + Wj,)/\/2 and vt = {vy — Vx)/\/2 
for the velocities along the diagonals. This gives 



9\\{rs) ^g{rs)+gyx{rs), 
9±{rs) ^g{rs)-gyx(rs), 



where we have also made use of the symmetry g^y = 
gyx which follows from considering the transformation 



— r, followed by a translation: Vx{Q)vy{Y) 



-> 



{-Vx{0))[-Vy{-Y)) = Vx{Y)Vy{Q) = Vy{0)Vx(v). 

Figure 6 shows the mixed correlations along the two 
diagonal directions. The thick arrows in the inset illus- 
trate the velocity fields one might expect as an effect of 
a particle with Vy > Q which is u^^ > for separations in 
the s direction and w^; < in the t direction. Since we ex- 
pect Uj; < along t, the correlations in that direction are 
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FIG. 6. The mixed correlation functions gyx{rs) and 
—gyx{rt). The difference between the mixed correlation func- 
tions (apart from the trivially different sign) is yet an example 
of the asymmetry in the sheared system. The arrows in the 
inset illustrate the rotational velocity fields which have the 
effect that a velocity Vy{0) > at the origin on the average 
gives Vx{rs) > and Vx{rt) < 0. 



shown with the opposite sign. Note that both functions 
start out at gyx{0) = {vxVy). The correlations then grow 
above this r = value and actually become stronger in 
the t direction (though with opposite sign) than along s. 
For a comparison, the dashed lines are the correlations 
of Vy along the same directions. This difference between 
5j/a;(^s) and —gyxift) is yet another example of a broken 
symmetry in the shearing system. 



IV. ROTATIONAL ASYMMETRY 

The previous section gave evidence for the importance 
of whirls in the velocity field. The inset of Fig. 6 shows 
two such whirls that are a consequence of Wy(0) > 0, and 
there will also be whirls with the opposite orientation 
that will contribute to the correlation functions in much 
the same way. The question we now like to address is 
whether they contribute equally much or not. Is the sys- 
tem symmetric when considering whirls with clockwise 
and counter-clockwise rotations, respectively? 

Since the shearing by itself introduces a rotational 
field, it could at first seem obvious that this symmetry is 
broken, but that is not correct. Our velocity correlations 
are calculated from the non-affine velocities and the net 
rotation is zero in the non-afRne velocity field. The sys- 
tem could therefore in principle well be symmetric with 
respect to these different directions of rotation. 



A. Asymmetric correlations 

There are at least two different ways to motivate the 
new correlation function that we are about to introduce. 
The first is to note that the direction of the rotational 
fields in the inset of Fig. 6 depend on the sign of Vy{0). 



With the opposite sign of Vy{0), Vx{ri) and Vx{rs) would 
also (typically) change sign and the rotations would be 
in the opposite directions. 

A second point of departure is to consider the fact that 
the correlation function defined in Eq. (1) is symmetric 
under the transformation x — > ~x whereas the system 
itself is only symmetric under the interchange of both 
X and y. This suggests that some information is lost 
when calculating the standard correlation function and, 
furthermore, that a guiding principle in the definition of 
an ideal correlation function is that it should have the 
same symmetry as the system. 

It is then possible to combine both these lines of 
thought and construct a function by restricting the av- 
erage in Eq. (1) to only include terms with Vy{r') > 0. 
For the transverse correlation function with r along the 
X direction, this becomes 



9a 



{vy{r')vy{r' + xx)) ^^1^^,-1^^. 



(2) 



This expression may further be generalized to including 
particles with both signs of Vy (r) by letting the direction 
of the separation {+xx or ~xx) depend on the sign of 
the velocity. 
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Vy{r')vy{r' +xx), Vy{r') > 0, 
Vy{r')vy{v' - xx), Vy{r') < 0. 



(3) 
The normalization A'term is the number of terms in the 
sum. Written this way it becomes clear that g^{x) in- 
deed has the desired symmetry properties. The ratio- 
nal for this new asymmetric function is also discussed 
in conjunction with Fig. 11 below. Note also that the 
symmetric function is related to g±{x) through 
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FIG. 7. Development of the asymmetry in G\(x). The func- 
tion is nearly symmetric at (^ = 0.82 but develops a very pro- 
nounced asymmetry as the density is increased above (/)j. The 
solid lines are the symmetric Gx. The shear rate is 7 = 10~ . 



gi.{x) = -^g\{x) +5+(-a;)] 



The quantities shown in the figures below (both symmet- 
ric and asymmetric) are, G{x) = g{x)/g{Q)^ normalized 
such that G(0) = 1. 

Figure 7 shows G\{x) at four different densities both 
above and below (/)j, again from simulations with 7 — 
10~^. For comparison, the symmetrized function G±{x) 
is given by the solid line. At </> = 0.82, panel (a), G'\_{x) 
is almost symmetric. The other panels show that the 
asymmetry grows with increasing ^, and &t (j) = 0.90 well 
above 0j, panel (d), the asymmetry is very pronounced. 
The a; < part shows a sharp dip at x sa — 18 whereas 
the a; > part has a rather shallow minimum at a; « 79. 

Another way to illustrate the growth of the asymmetry 
is through the position of the minima. Fig. 8 shows (.^ 
and (.- , the absolute value of the position of the minima 
for x > and a; < 0, respectively, together with tgym 
from the minimum of the symmetrized function. The 
asymmetry grows rapidly above = 0.84 and the behav- 
ior of £-|_ and £_ turn out to be very different at higher 
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FIG. 8. Location of the respective minima of the correlation 
functions. £+ and I- are the positions of the minima of Gt 
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FIG. 9. The correlation function G+{x) at = 0.8433 « 0j 
and several different shear rates down to 7 = 10~*. The 
asymmetry increases slowly with decreasing shear rate. 



Fig. 9 is the same quantity obtained at (j) — 0.8433 « 
4>j [13] with different shear rates. It is here found that 
the asymmetry grows with decreasing shear rate. The 
algebraic increase of £+, which suggest a divergence in 
the hmit of vanishing shear rate, is clear from Fig. 10 
whereas both £_ and £sym grow less rapidly. 



B. Origin of the rotational asymmetry 

We now turn to the question of the origin of the asym- 
metry in g]^{x) and will argue that it is linked to the 
plastic processes. The reason for this is the manifest 
asymmetry of the velocity profile of an elementary plas- 
tic event, as shown in the inset of Fig. 11 [9]. Note that 
the orientation of this velocity profile with its quadrupo- 
lar structure is dictated by the direction of the shear. 
The velocity field corresponds to a compression along the 
±{x — y) direction together with an expansion in the or- 
thogonal ±{x + y) directions, which is equivalent to a 
simple shear. This implies that the mirror version of 
such an event is expected to be much less common, if at 
all present, and it is this effect that causes the asymme- 
try. Note that this agrees well with the enhanced anti- 




FIG. 10. Position of the minima oi g'\^{x) and g±, respectively, 
a,t <j) — 0.8433 ~ (j)j and four different shear rates. Note 
the algebraic increase oi £+ ~ 7" ■ ^ which would imply a 
divergence at vanishing 7. 



correlation at a; < 0. When a particle has Vy > due 
to a plastic event one would expect that there should be 
one or more other particles with Vy < 0, and the inset 
of Fig. 11 shows that they should be found in the —x 
direction. 

To further check this idea we have tried to separate the 
contribution to the correlation function from the plastic 
events from the rest, i.e. we have calculated gj^ sepa- 
rately for increasing and decreasing total energy [6], re- 
spectively. That study, which was done for N — 4096 
particles, did indeed give evidence (not shown) that the 
part related to a decrease in total energy was the more 
asymmetric one. However, since the total energy is a 
global quantity that kind of analysis isn't a very sensi- 
tive one, at least not at finite shear rates. In a large sys- 
tem one would expect some regions to be characterized 
by plastic events and a local decrease in energy whereas 
the motion in other regions is elastic, but we still have to 
classify the whole system as either plastic or elastic. 

Instead of splitting up the contributions based on the 
change in the total energy we now consider the change in 
the local energy. That is done by classifying each term 
in Eq. (3) as "fast" or "slow", corresponding to plas- 
tic and elastic respectively. Since the change in total 
energy is AE = L'^a'j — X^i '^? [6]; the change in the 
local energy is related to the average v^ in a certain re- 
gion. On the average, one expects the energy to decrease 
locally if v^ > (v^). As a reasonable and more sensi- 
tive way to split gj_{x) into two different contributions, 
dl = .9±°" + 5'l''* (we drop the "4-" in the "slow" and 
"fast" terms to simplify the notation) we therefore clas- 
sify each of the terms in Eq. (3) according to the mag- 
nitude of the velocities. If both v(r') and v(r' -I- xx) are 
"slow", (i.e. they both obey v^ < (v^)), the term con- 
tributes to g'f"^ whereas it contributes to g^''^' otherwise, 
i.e. if at least one of the particles is "fast" . 

Figure 11 shows the splitting of the total G~^{x) (open 
circles) into slow and fast parts, respectively. Note that 




symmetries due to the shearing. We first find that 
{vxVy) y^ for individual particles and that this cor- 
relation depends strongly on density. At low densities 
the particle motion is preferably along the force chains 
whereas it is preferably perpendicular to the force chains 
around (j)j. We then examine how the total velocity 
correlation g{r) depends on the direction of r. Rather 
surprisingly the correlation along the two diagonals, cor- 
responding to the direction of compression and dilation 
respectively, are almost identical. The decay along the 
diagonals is non-monotonic, in contrast to the monotonic 
decay along the main (x and y) directions. 



FIG. 11. Splitting of the correlation function into two dif- 
ferent contributions at 4> = 0.8433 and 7 — 10^^. The total 
G^{x) (open circles) is split into contributions from slow and 
fast particles, respectively. It is clear that the asymmetry is 
largely an effect of the fast particles which we relate to the 
plastic processes. The inset shows an idealized velocity profile 
for a typical plastic event with quadrupolar structure. Note 
that this velocity field is a consequence of our shearing geom- 
etry which is equivalent to a compression along the ±{x — y) 
directions together with an dilation in the orthogonal ±(x + y) 
directions. 



the contribution from the slow particles (solid dots) is 
almost symmetric whereas there is a very pronounced 
asymmetry from the fast particles (open squares). This 
is therefore evidence that the asymmetry is caused by the 
fast particles and that this typically is related to a local 
drop in energy which is often related to a plastic event. 



We then argue that the usual correlation functions 
are more symmetric than the system itself and define 
a less symmetric velocity correlation function that also 
may be used to probe the differences between clockwise 
and counter-clockwise rotations. This function is asym- 
metric with respect to x, and this is an asymmetry that 
increases rather dramatically when either the density in- 
creases above (pj or the shear rate decreases at fixed 
(j) ■^ (f)j. We attribute the asymmetry to elementary plas- 
tic events with quadrupolar symmetry and their orienta- 
tion dictated by the direction of the shear as shown in 
the inset of Fig. 11. 



V. SUMMARY 

To summarize, we have examined velocity correlations 
in a sheared system with emphasis on the breaking of 



I thank S. Teitel for helpful discussions and critical 
reading of an earlier version of the manuscript. This work 
was supported by the Swedish Research Council and the 
Swedish National Infrastructure for Computing. 



[1] A. J. Liu and S. R. Nagel, Nature (London) 396, 21 [7] 

(1998). 
[2] P. Olsson and S. Teitel, Phys. Rev. Lett. 99, 178001 [8] 

(2007). 
[3] T. Hatano, J. Phys. Soc. Jpn. 77, 123002 (2008). [9] 

[4] M. Otsuki and H. Hayakawa, Phys. Rev. E 80, 011308 

(Jul 2009). [10] 

[5] B. P. Tighe, E. Woldhuis, J. J. Remmers, W. van 

Saarloos, and M. van Hecke, "Model for the scaling of [11] 

stresses and fluctuations in flows near jamming," (2010), 

arXiv: 1003. 1268. [12] 

[6] I. K. Ono, S. Tewari, S. A. Langer, and A. J. Liu, Phys. [13] 

Rev. E 67, 061503 (2003). 



G. Lois, A. Lemaitre, and J. M. Carlson, Phys. Rev. E 

76, 021302 (2007). 

T. S. Majmudar and R. P. Behringer, Nature 453, 1079 

(2005). 

C. Maloney and A. Lemaitre, Phys. Rev. Lett. 93, 016001 

(2004). 

C. S. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, 
Phys. Rev. E 68, 011306 (2003). 

D. J. Evans and G. P. Morriss, Statistical Mechanics of 
Nonequilibnum Liquids (Academic Press, London, 1990). 
D. J. Durian, Phys. Rev. Lett. 75, 4780 (Dec 1995). 

C. Heussinger and J.-L. Barrat, Phys. Rev. Lett. 102, 
218303 (May 2009). 



